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Radiation-Driven Warping. II. Non-Isothermal Disks 

Philip R. Maloney 1 , Mitchell C. Begelman 2,3 and Michael A. Nowak 2 

ABSTRACT 

Recent work by Pringle and by Maloney, Begelman & Pringle has shown that 
geometrically thin, optically thick, accretion disks are unstable to warping driven 
by radiation torque from the central source. This work was confined to isothermal 
(i.e., surface density S oc R~ 3 / 2 ) disks. In this paper we generalize the study of 
radiation-driven warping to include general power-law surface density distributions, 
£ oc R~ s . We consider the range 5 = 3/2 (the isothermal case) to 5 = —3/2, which 
corresponds to a radiation-pressure-supported disk; this spans the range of surface 
density distributions likely to be found in real astrophysical disks. In all cases there 
are an infinite number of zero-crossing solutions (i.e., solutions that cross the equator), 
which are the physically relevant modes if the outer boundary of the disk is required 
to lie in a specified plane. However, unlike the isothermal disk, which is the degenerate 
case, the frequency eigenvalues for 5 3/2 are all distinct. In all cases the location of 
the zero moves outward from the steady-state (pure precession) value with increasing 
growth rate; thus there is a critical minimum size for unstable disks. Modes with zeros 
at smaller radii are damped. The critical radius and the steady-state precession rate 
depend only weakly on 5. 

An additional analytic solution has been found for 6 = 1. The case 5 = 1 divides 
the solutions into two qualitatively different regimes. For 6 > 1, the fastest-growing 
modes have maximum warp amplitude, /3 max , close to the disk outer edge, and the 
ratio of /3 max to the warp amplitude at the disk inner edge, /3 G , is 3> 1. For 5 < 1, 
/3max//3o — 1) and the warp maximum steadily approaches the origin as 5 decreases. 
This implies that nonlinear effects must be important if the warp extends to the disk 
inner edge for S > 1, but for 5 < 1 nonlinearity will be important only if the warp 
amplitude is large at the origin. Because of this qualitative difference in the shapes of 
the warps, the effects of shadowing of the central source by the warp will also be very 
different in the two regimes of 5. This has important implications for radiation-driven 
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warping in X-ray binaries, for which the value of 5 characterizing the disk is likely to 
be less than unity. 

In real accretion disks the outer boundary condition is likely to be different from 
the zero-crossing condition that we have assumed. In accretion disks around massive 
black holes in active galactic nuclei, the disk will probably become optically thin before 
the outer disk boundary is reached, while in X-ray binaries, there will be an outer 
disk region (outside the circularization radius) in which the inflow velocity is zero but 
angular momentum is still transported. We show that in both these cases the solutions 
are similar to the zero-crossing eigenfunctions. 



Subject headings: accretion disks - instabilities - galaxies: individual (NGC 4258) - 
stars: individual (Her X-l, SS 433) - X-rays: stars 



1. Introduction 

Evidence for warped, precessing accretion disks in astrophysical systems ranging from X-ray 
binaries to active galactic nuclei has steadily accumulated over the last two decades (see Maloney 
& Begelman 1997a, and references therein). The origin and maintainence of such warped disks 
has until recently stood as an unsolved theoretical problem. While it is possible, for example, to 
generate non-planar modes with m = 1 symmetry in thin, relativistic disks flKato 1990| ; |Kato fc 



Honma 1991), these modes only exist at small radii (R 10 Schwarzschild radii), since they rely on 
trapping of the modes in the non-Newtonian region of the potential. However, an important clue 
was provided by |Petterson (1977) , who pointed out that in an optically thick disk with a central 



source of luminosity, the pressure resulting from re-radiation of the intercepted flux will produce a 
net torque if the disk is warped. Almost twenty years were to pass before it was recognized that 
radiation pressure torque actually leads to a warping instability. Pringle (1996)| (P96) showed that, 



for the special case in which the disk surface density E oc R~ 3 / 2 (corresponding to an isothermal 
disk in the usual a— disk formalism, with disk viscosity v oc R 3 / 2 ), even an initially planar disk is 
unstable to warping by this mechanism. Pringle solved the linearized twisted disk equations in 
this case using a WKB approximation. Maloney, Begelman & Pringle (1996, Paper I, hereafter 
MBP) found exact solutions to the linearized twist equations, and demonstrated the importance 
of the outer boundary condition for determining the growth rates of the unstable modes. 

These previous works all specialized to the case of an isothermal disk, which simplifies the 
twist equations. While this may be a reasonable approximation for some astrophysical disks 
(e.g., the masing molecular disk in NGC 4258; see MBP), there are many other systems, such 
as accretion disks in X-ray binary systems, where this is likely to be a poor assumption. In 
this paper we extend the work of MBP by considering disks with power-law surface density 
profiles, E oc R~ s . We consider the range —3/2 < 5 < 3/2: the lower limit corresponds to a 



radiation-pressure-supported disk (e.g., Frank, King & Raine 1992), while the upper limit is 
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the isothermal value (MBP). Within the limitations of assuming a constant power-law for the 
surface density, this spans the probable range of surface density laws relevant to real astrophysical 
accretion disks. For example, the standard Shakura-Sunyaev gas pressure-supported disk is 
characterized by 5 = 0.75 ( |5hakura &: Sunyaev 1973| ). 

In §2 we discuss the twist equation, including the effect of radiation torque, and cast it 
into a more convenient form. We solve the equation numerically in §3 and discuss both the 
time-dependent and steady-state solutions. As in the isothermal case, the outer boundary 
condition is crucial for determining the stability of the disk and the growth rates of the unstable 
modes. In §4 we discuss the important issue of the appropriate outer boundary condition for 
accretion disks around stellar-mass objects and AGN. Finally, in §5 we discuss the implications of 
the results and their application to real accretion disks. 



2. The Disk Evolution Equation 

As in P96 and MBP, we adopt a Cartesian coordinate system with the Z— axis aligned with the 
normal to the unwarped disk, and define (3 to be the local angle of tilt of the disk axis with respect 
to the Z— axis, while 7 — ir/2 is the angle between the descending line of nodes and the A— axis. 
The equation governing the evolution of the local tilt vector, l(R, t) = (sin (3 cos 7, sin (3 sin 7, cos (3) , 
including the radiation torque term, is then (P96) 
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where V\ and V2 are the disk viscosity in the azimuthal and vertical directions (with a ratio 
i] = vilv\ that is assumed constant but not necessarily unity), primes denote derivatives with 
respect to R, £ is the disk surface density, f2 the Keplerian angular velocity, Vr the radial inflow 
speed, and the radiation torque term 

where dG is the torque exerted on a ring of width dR and radius R and is given by equation (2.18) 
of P96. 

As in MBP, we assume the disk viscosity can be written is\(R) = v {R/R ) 5 , where R is 
an arbitrary fiducial radius, but we now allow the radial power-law to be arbitrary, rather than 
specializing to the case 5 = 3/2. In a steady-state disk far from the boundaries, S = M/2>ttv\^ 
so this radial dependence of the viscosity translates directly into a power-law surface density, 
£ oc R~ s . Furthermore, in a steady disk the radial inflow velocity Vr = v±Q' /Q, and so the first 
two terms inside the brackets on the LHS of equation (1) cancel, while the third term becomes 

{ER3n) '-f 3 *U. (3) 



T,R 3 n \2 J R 
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We then linearize to obtain the more general version of equation (1) of MBP for W = (3e n , 

dW /3u 2 „\ dW (l d 2 W\ 

' ' ^UrtT > ( 4 ) 



dt \AR ) OR \2 OR 2 J ' 
where the radiation torque term is 

r = n^ERXte ^ 
with L the luminosity of the central source. Assuming an accretion-fueled source with radiative 
efficiency e = L/Mc 2 , this term can be written as 

^R S - 1/2 = ToR 5 - 1 / 2 (6) 



where R s is the Schwarzschild radius (cf. equation [3] of MBP). Transforming to the new radius 
variable x = aR 1 ^ 2 , using equation (6) for T, and Fourier transforming with respect to time, the 
twist evolution equation becomes 



dt 8R 5 n 



d 2 W ( iV2e \ dW 
X dx2 V vaRl^J dx 



(7) 



where in general a is complex. We set a = V2e/R^ 2 r] (note that this differs from the definition 
of radius variable z in MBP only in making x real, rather than pure imaginary). Since R Q is 
arbitrary, we now define R Q so that x(R Q ) = 1; v is thus the value of v\ at x = 1. Therefore 

Ro = a' 2 = M (8) 

and the coefficient of the x 2S ~ 3 term on the righthand side of equation (7) becomes u e 4 /2R 2 rj 3 . 
Finally, we nondimensionalize the eigenfrequency a by defining 

to obtain the final form of the twist equation: 

d 2 W dW , or 

x^r- + (2 - ix)^- - iax 3 ~ 2S W = . (10) 
ox 1 ox 

For 5 = 3/2, equation (10) reduces to Kummer's equation, as discussed in MBP. (Note that 
the coefficient of W in equation [9] of MBP, 2a/T , is identical to a as defined here.) Equation 
(10) can thus be regarded as a modified Kummer's equation. Unfortunately, this equation is in 
general analytically intractable, so that numerical solution is necessary. We discuss the asymptotic 
behavior of the solutions in Appendix A. For the physically relevant solutions (those which exhibit 
zero-crossings), there are in all cases an infinite number of zero-crossing eigenfunctions. As we will 
see below, the special case 5 = 3/2 is degenerate: the real parts of the eigenvalues are all identical, 
with Re(cr) = 1. For all the other values of S that we consider the eigenvalues are distinct, and 
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each eigenfunction has at most one zero. In addition, the real part of a, 5y, must be greater than 
zero, which means that the direction of precession of the warp must be the same as the direction 
of disk rotation (i.e., prograde: see Appendix B). This is a simple consequence of the overall 
shape of the warp: it is easy to show that for a disk in which the gradient in the tilt /3' > 0, the 
direction of precession due to the radiation torque is retrograde, while if j3' < the precession will 
be prograde. Since the solutions are constrained to return to the plane at the outer boundary, 
is always negative at large radius, and this dominates the direction of the induced precession. 

We have found one additional analytic solution to equation (10), for 5 = 1 (Appendix C). We 
discuss this and the numerical solutions to equation (10) for both steady-state and unstable modes 
in the following section. 

3. Solutions of the Twist Equation 

Real astrophysical accretion disks will generally be unwarped beyond some maximum radius 
i? max , either because the disk becomes optically thin to the incident or re-emitted radiation (see 
the discussion in MBP), so that the disk must eventually return to the initial plane, or because 
this is forced in some other manner by the physical outer boundary condition (e.g., if the accretion 
disk is fed by material from a companion star; see §4.2). We therefore consider solutions of the 
twist equation for which the disk returns to the original plane Z = at some radius. This outer 
boundary condition is not rigorously correct, as we discuss in §4, but it is in general an excellent 
approximation to the true outer boundary condition, and the zero-crossing eigenfunctions are very 
useful for understanding the behavior of the warping instability. 

Equation (10) is a second-order differential equation, and therefore requires two boundary 
conditions to specify the solution. In principle it could be solved as a two-point boundary value 
problem, except that the location of the outer boundary is not necessarily known a priori. It is 
computationally most convenient to solve it as an initial-value problem, making an initial guess 
for the value of a and then iterating to find the solution that goes to zero. We separate equation 
(10) into a coupled pair of equations for the real and imaginary parts of W, and integrate outward 
from the origin using a fifth-order Runge-Kutta scheme ( press et al. 1992] ). 

We require solutions that are regular at the origin. As x —* 0, the leading terms of equation 
(10) are 

xW" + 2W' = ix 3 ~ 25 aW (11) 

where primes now denote derivatives with respect to x. Adopting the boundary condition 
W(0) = 1, the leading behavior of (11) as x — > gives 

W'^^-.x 3 - 25 as x^0. (12) 
5 — 2o 

Furthermore, we require that zero torque be exerted on the disk at the origin, which requires that 
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x 2 W — > as x — > 0. From equation (12), it is evident that this second boundary condition is 
satisfied for any 5 < 5/2. 

Numerical solutions to the twist equation are presented below. However, we first note that 
there are simple scaling relations that are generic features of the instability: 

(1) The radiation torque per unit area is V ~ (L/4itR 2 c) x R, while the angular momentum 
density is I ~ QR 2 T, ~ QR 2 M /3irvi, where £1 is the angular velocity at radius R in the disk, 
M is the mass accretion rate, and v\ is the usual kinematic viscosity. The radiation torque 
timescale is thus given by 

OR 3 Mc 2 FIR 3 , , , 

trad r~ = e" 1 13 

which, for an accretion- fueled source, depends only on the accretion efficiency e = L/Mc 2 
and not on the luminosity L and the mass accretion rate individually (MBP; Maloney & 
Begelman 1997a). 

(2) The viscous timescale is t v - lsc ~ R/Vr = 2R 2 /3ui, so the ratio of viscous to radiation torque 
timescales is given by 

which is independent of the form of the viscosity law and depends only on e and the radius 
normalized to the gravitational radius. Thus the radiation torque always wins out at large 
radii, i.e., for R > e~ 2 R g , but viscosity always dominates near the center. This is why 
the disk will be flat (but in general will have non-zero tilt) at small radii. Note that this 
equation also implies that t v iscArad ~ 1 at x = 1. 

(3) As is apparent from equations (7) - (9), altering the values of r\ and e while keeping the ratio 
rj/e fixed will affect the growth rates, but not the shapes of the solutions, i.e., the value of x 
is unaffected. 




3.1. Steady-State Solutions 

We first discuss the steady-state solutions to the twist equations, for which <jj = Im(cr) = 0. 

Figure 1 shows the location of the zero, x Q , for the first ten zero-crossing eigenfunctions, as a 
function of surface density index 5. As noted earlier, for 5 ^ 3/2 the eigenfunction corresponding 
to each eigenvalue has only a single zero (Appendix B). We use "order" to specify how many 
eigenfunctions have their zero at values of x smaller than or equal to the eigenfunction in question; 
thus the eigenfunction with the smallest value of x Q is the first-order eigenfunction, that with the 
next smallest is the second-order eigenfunction, and so forth. The behavior of the higher-order 
eigenfunctions is remarkably complex. 
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In Figure 2 we plot the normalized real eigenvalues, ay, for the first ten order eigenfunctions. 
The first order eigenfunction has the largest eigenvalue for all values of 5. The degeneracy of the 
case 5 = 3/2 discussed in MBP is immediately apparent, with ay = 1 for all the eigenfunctions. 
For 5 / 3/2, the degeneracy is lifted, with increasing separation of the eigenvalues as 5 decreases. 
Closer examination shows that the behavior becomes quite complex for the higher-order zeros 
with decreasing 5, as merging of eigenvalues occurs. 

The location of the first-order zero as a function of S has an important physical significance. 
In real astrophysical disks, shadowing effects probably make the higher-order eigenfunctions 
unimportant. Just as for the isothermal case, for the growing modes the location of the zero moves 
outward from the steady-state value as the growth rate increases (see below); all the modes with 
zeros at smaller radius are damped. Thus the location of the first zero marks the critical boundary 
for disk stability. Accretion disks larger than 



are unstable to warping by radiation pressure, where the critical value x cr is equal to 2tt for 
5 = 3/2 and x cr — 4.8917T for 5 = —3/2. Figure 1 shows that x cr increases smoothly as 5 decreases. 
Disks with outer boundaries smaller than R cr are stable against radiation-driven warping. From 
equation (15), this critical radius scales as e™ 2 , so that accretion disks in low-efficiency systems 
will not be unstable unless the disks are extremely large. (Furthermore, the precession and growth 
timescales for the instabilty will be extremely long if e is very small; cf. the discussion of equations 
[30] and [31] in §5.) 

In Figure 3 we plot the tilt [3 as a function of radial variable x for the first-order steady-state 
solutions, for several values of 5. For clarity they have been plotted only to the zero of the mode. 
In all cases the maximum tilt is at the origin (note that /3 max is arbitrary and has been taken to 
be unity); as 5 decreases from 3/2 the zero moves to larger x and an increasing fraction of the disk 

has (3 ~ (3 max . 

The steady-state modes have shapes that are time-independent in a frame rotating with 
angular velocity a r ; physically, these are purely precessing modes. The frequency a r is related 
to the dimensionless frequency ay by equation (9). As can be seen from Figure 2, ay for the 
first-order zero drops by six orders of magnitude from 5 = 3/2 to 5 = —3/2. However, to put the 
results in physical terms, note that from the definition of x, we can rewrite equation (9) as 



Defining the viscous timescale as before as 




(15) 



^visc(-R) — „ j 

we can then express a in terms of the viscous timescale at the critical radius, 



(17) 



a = 




(18) 
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where x cr is the corresponding value of x (cf. equation [15]). For S = —3/2, x^~ 2S a r = 176.8, 
while for 5 = +3/2 this quantity equals 2ir, so that, in terms of the viscous timescale at R cr , the 
variation in o r is only a factor of ~ 28 over the entire range of 6. Alternatively, a can be expressed 
in terms of the viscous timescale at the outer edge of the disk (i.e., the location of the zero), with 
x cr and R cr replaced by x and R(x ). (In fact, the variation in a r is smaller than this, because 
for a fixed value of the viscosity, the increase in R CI with decreasing 5 causes t v i sc (R cr ) to increase, 
offsetting the variation in xf~ 25 a r ; see §5.) 

3.2. Unstable Solutions 

We now consider solutions that evolve with time, i.e., &i / 0. As for the isothermal disks 
analyzed in MBP, the only modes that have zeros at radii smaller than x cr , as discussed above, 
are damped, with a, > 0, and we do not consider them further. 

For the case of an isothermal disk considered by MBP, the zeros move steadily outward 
with increasing growth rate (i.e., increasing — <7j), and there is no limit to the growth rate. For 
the non-degenerate cases considered here, the zeros also move outward with increasing — <7j, as 
expected from the scalings discussed in §3. However, since the degeneracy of the eigenvalues has 
been lifted, there is now a maximum growth rate for a given order mode. In particular, in Figure 
4, we plot the location of the first-order zero, x a , (as discussed above, this is probably the only 
physically relevant eigenfunction) as a function of —5*,, for 5 = —3/2 to 1.45, in steps of 0.05 in S. 
The systematic behavior of the location of the zero is quite striking, and there is a clear transition 
across 5 = 1. The maximum possible (normalized) growth rate — d{ for the first-order zero steadily 
decreases as 5 decreases. But the location of the zero at this maximum growth rate, — oy max , 
moves to larger x as 5 decreases from 1.45 to 1.00, even though — oy max systematically declines 
with decreasing 5. 

The case 5 = 1 is special. There is an analytic solution in this discussed in Appendix 

C. The maximum growth rate is — oy m ax = 0.25; at this value of <7j, the zero has moved to x = oo. 
(The curve plotted in Figure 4 ends at — <5y max = 0.249.) This value of 5 is also special in that 
all of the eigenvalues merge as ai — ► ay max , so that the zeros of all orders move to infinity at 
— <5"i,max = 0.25. Figure 5a plots the location of the first three order zeros for 8 = 1 as a function of 
— <7i, while Figure 5b shows the real parts of the corresponding eigenvalues. 

For 5 < 1, the maximum x— value for the first-order zero decreases rapidly as 5 decreases, 
reaching a minimum at 5 ~ 0.15, and then slowly increasing as 6 approaches —3/2 (see Figure 4). 
The ratio x (ai^ mSLX )/x (ai = 0) decreases as 5 decreases, so that the range of instability (out to 
the first zero) decreases with 5 in this range. 

In Figure 6 we plot the real part of the eigenvalue, ay, versus the imaginary part for the 
growing modes, for the same values of 5 as in Figure 4. In all cases ay shows a characteristic 
decrease as oy max is approached. As in Figure 4, the curve for S = 1 is plotted only to — ai = 0.249; 



-9- 



a r — > as — <j\ — > 0.25. The magnitude of the drop in <5y declines as <5 approaches 3/2; for the 
isothermal 1, independent of <7j (MBP). 

There is a dramatic change in the behavior of the unstable solutions across 5 = 1. In Figure 
7 we have plotted the shape of the eigenfunctions, i.e., (3 as a function of x, for several values of 
5 > 1, for the fastest-growing modes, with Oi = Oj jmax . The eigenfunctions have been normalized 
so that the maximum value of (3 is unity; the normalization factor (i.e., l//3 max ) is very small, of 
order 10 -6 . The behavior of the eigenfunctions is very similar to the fast-growing modes of the 
isothermal disks discussed in MBP: the warp has its maximum close to the outer boundary of the 
disk (assumed to be coincident with the zero), and the amplitude of the warp at this maximum is 
very large compared to that at x = (approximately 10 6 for the modes shown in Figure 7). The 
shape of the fastest-growing modes for 5 = 1 are qualitatively similar to the 5 > 1 modes; the only 
difference is that x Q approaches oo (rather than a finite value) as a approaches <7j imax . 

In Figure 8 we similarly plot the normalized eigenfunctions for the fastest-growing modes, 
for several values of 5 < 1. The change in the nature of the eigenfunctions is marked: even for 
5 = 0.95, the amplitude of (3 at the maximum is only ~ 2.5 times the value of (3 at x = 0, (3 , 
and as 5 decreases /3 max /A> declines toward unity. For 5 ^ 0.5, the eigenfunctions converge to 
an essentially constant shape, with a plateau of constant (3 extending from the origin to x ~ 5, 
followed by a decrease to the zero. This change in the behavior of the eigenfunctions across 
5 = 1 has important implications, as we discuss below. Note also that this applies only to the 
unstable modes; the steady-state solutions do not exhibit any change in behavior across 5 = 1, as 
is apparent from Figure 3. 

4. Outer Boundary Conditions 

In section 3, we solved equation (10) with the requirement that the solution go to zero at 
the disk outer boundary, which is at some finite radius. While this produces a well-defined set 
of solutions whose behavior is very useful for understanding the behavior of the instability, in 
real astrophysical disks the true outer boundary condition is likely to be different. Furthermore, 
the correct outer boundary condition for accretion disks around stellar-mass compact objects 
(neutron stars and black holes) in X-ray binaries will differ from that appropriate for accretion 
disks around massive black holes in active galactic nuclei. In this section we discuss the choice 
of outer boundary condition and how the solutions in these cases differ from the zero-crossing 
solutions discussed above. 



4.1. Active Galactic Nuclei 

It is expected that accretion disks will generally satisfy the requirement of optical thickness, 
and will therefore be subject to the radiation- warping instability (see §5 and Appendix D). 
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However, accretion disks surrounding massive black holes in active galactic nuclei are likely to 
become optically thin (to the re-emitted radiation) before the physical boundary of the accretion 
disk is reached. (This is true for 5 > 0, so that the disk surface density decreases with increasing 
radius; this is undoubtedly the case at the relevant radii of AGN accretion disks.) Once the disk 
becomes optically thin, the instability ceases to operate, as the re-emitted radiation no longer 
exerts a torque on the disk (cf. MBP). Therefore, equation (10) does not describe the behavior 
of the disk beyond the optically thin radius .Rthin- However, viscosity continues to operate, so 
that twist angular momentum generated by the instability can be transported beyond -Rthin- We 
therefore have to solve a pair of equations, with the solutions matched at -Rthin : interior to -Rthin > 
the equation governing the dynamics of the disk is equation (10), while exterior to -Rthin the disk 
is described by equation (10) minus the torque term —ixW. Both W and W' must be continuous 
at -Rthin! the outer disk solution must also obey W — > as R — > oo. In a real accretion disk, the 
radiation torque will presumably decline smoothly to zero rather than switching off abruptly, but 
as long as this occurs in a radial thickness AR <C R we expect this to have little effect on the 
solutions. 

For 6 = 3/2 and 5 = 1 the outer disk equation can be solved analytically; the solutions are 
modified Bessel functions {K\ and K 1 / 2 , respectively). In general it must be solved numerically. 
For a specified value of -Rthin we solve both the inner and outer disk equations and iterate to find 
the solution that satisfies the matching condition at -Rthin- (In practical terms, the constraint that 
must be satisfied is that (3 1 be continuous at .Rthin, since (3 can be scaled arbitrarily and 7 does 
not enter into equation (10), only 7' and 7"; the value of 7 can always be matched by an arbitrary 
rotation of the outer disk solution.) Since the parameter space to be considered is very large, we 
consider only a few representative solutions for 5 = 3/2, 1.25, and 0.75; the behavior of solutions 
for 5 < 0.75 is very similar to that of the 5 = 0.75 solutions. 

We first consider the modifications to the steady-state, first-order solutions shown in Figure 3. 
We take -Rthin to be the zero of these modes, and then find the solutions that satisfy the boundary 
conditions at .Rthin- In Figure 9 we plot the real versus the imaginary parts of the eigenvalues for 
these solutions. Not surprisingly, there is now a continuum of solutions rather than a single mode, 
as it is easier to satisfy the boundary condition (3 — > as x — ► 00 rather than at a finite value of x. 
The solid lines correspond to growing modes (cr, < 0) while the dotted portions of the curves are 
damped. The shapes of the solution curves in the (ay, o"j)-plane are very similar for all values of 
5, with a pair ay-values for a given value of o^. In all cases there is a well-defined maximum (and 
minimum) growth rate at which the two values of ay merge; the precession rate ay corresponding 
to this fastest-growing mode differs from that of the zero-crossing solution by ~ 2%, 8.5%, and 
13% for 5 = 1.5, 1.25, and 0.75, respectively. 

In Figure 10 we plot the tilt (3, normalized to a maximum of unity as before, for the 
steady-state zero-crossing modes (dotted lines) and the corresponding fastest-growing modes from 
Figure 9, with -Rthin = x o (solid lines). The behavior of the "AGN" modes is what one would 
expect: since the growth rates are non-zero, the maximum of the tilt has moved outward from the 



- 11 - 



origin, and the solutions decline smoothly to zero, with W' — > as x — ► oo. Since the differences 
in shape between the steady-state mode and growing modes are minor for 5 0.75, as discussed 
in §3, the differences between the zero-crossing modes and the optically-thin-boundary modes are 
minimal for the 5 = 0.75 case. We also note that beyond -Rthin j' becomes negative (i.e., the line 
of nodes describes a retrograde spiral), since in the absence of radiation torque viscosity simply 
causes the line of nodes to unwind. 

We now consider modifications to the rapidly growing modes. Figure 11 plots the precession 
rate versus the growth rate for rapidly growing modes for the three values of 5. For 5 = 1.25 and 
0.75 the zero-crossing modes used to set the value of i?thin have growth rates close to the maximum 
for first-order modes. For 5 = 1.50 there is no maximum growth rate (due to the degeneracy) 
and so the solution shown is simply a rapidly growing mode, with &i = —5. The behavior of the 
continuum of modes that satisfies the boundary conditions is very similar to those shown in Figure 
9, except that there are no damped modes in this case. Comparison with Figure 9 also shows 
that for 5 ^ 1.50, the width of the curve in the (ay, <7j)-plane is significantly smaller for these 
fast-growing modes. The precession frequencies of the maximally-growing modes are even closer 
to those of the zero-crossing modes than for those shown in Figure 9, differing by ~ 0.05%, 1%, 
and 8% in order of decreasing 5. 

In Figure 12 we plot (5 as a function of x, as in Figure 10, with the fastest-growing modes 
shown as solid lines and the corresponding zero-crossing modes shown as dashed lines. The 
differences between the zero-crossing solutions and the x(-Rthin) = x Q modes is even smaller than 
in Figure 10; the peak in (5(x) is displaced slightly outward and (3 declines smoothly to zero with 
increasing x. For the 5 > 1 modes, which peak well away from the origin, the width of the peak in 
(3(x) at half-maximum is negligibly different in the two cases. For 5 ?S 0.75, it is evident that the 
difference between the zero-crossing modes and the optically-thin-boundary modes is negligible for 
all growth rates. 



4.2. X-Ray Binaries 

In X-ray binary systems the situation is rather different. It is likely that the accretion 
disk remains optically thick to the physical outer boundary of the disk. However, the disk will 
again not obey equation (10) throughout, because the assumption of a power-law surface density 
cannot hold for the entire disk. What we have taken to be the outer boundary of the disk in 
the zero-crossing solutions corresponds to the circularization radius, .R c i rc , in a real X-ray binary 
system. The disk will actually extend to the tidal truncation radius R on t, which is larger by a 
factor of 2 — 3 for reasonable mass ratios (e.g., Frank, King fc Raine 1992] ). In this outer disk 



region, i? c irc < R < Rout, the radial velocity Vr = 0, since there is no flux of mass through this 
region of the disk, only angular momentum. Conservation of angular momentum then requires 
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that 

^=(^)^(%*) 1/2 . (19) 

To determine the surface density distribution in the outer disk thus requires an assumption about 
the viscosity. We assume that v\ = u out is constant in the outer disk; in fact, since i? ou t is only 
~ 2 — 3i? c i rc , the results are entirely insensitive to this assumption unless v is an extremely strong 
function of radius in the outer disk. The term in brackets on the LHS of equation (1) is then 



n' i (ER 3 n)' 

Vr - vi— - -i/ 2 - 



^(3-^7). (20) 



Q 2 ER 3 Q 

Fourier transforming with respect to time and linearizing as before, the twist equation becomes 

"77-3 2ir out 



r)W" + 



R ^out 



W-™?- = (21) 

^out 



where r out , defined as in equation (5), is constant since £i? 2 f2 is also constant under the 
assumption that v is independent of radius in the outer disk. This does not include the torque 
due to the companion star, which can be of vital importance for the disk modes in X-ray binary 
systems: as shown in Maloney & Begelman (1997b), the companion torque allows retrograde 
modes as well as prograde modes to exist. Since the orbital period is short compared to the 
viscous timescale in the disk, we average over azimuth and keep only the leading (quadrupole) 
term in the companion torque. This contributes a term —itu (R/R ) 3 ^ 2 W to the lefthand side of 
equation (21), where oj is the quadrupole precession frequency at the fiducial radius R a (Maloney 
& Begelman 1997b). 

By continuity, the viscosity i/ out and the radiation torque parameter r out in the outer disk are 
given by 

"out = ^l(-Rcirc) = V Q {Rdrc/Ro) (22) 

Tout = To-Rdr^ 2 • (23) 

Transforming to radius variable x and nondimensionalizing both a and uj by 2r/ 3 i?^/f e 4 (cf. 
equation [9]) the linearized twist equation for the outer disk becomes 

"2(77-3) i - 2 - 3 



xW" + 



- 1 



V 



W>-^W>=™(o + u, x*)w. (24) 

x circ X circ 



For typical neutron star X-ray binary parameters, x c i rc w 30 (assuming e w 0.1), while 
%out — \/3x c i rc . 

The value of W at the outer boundary is arbitrary, but we require that W = at i? ut; this 
implies that, in the absence of tidally-induced precession, there is no torque acting on the outer 
disk boundary. As in the case of the AGN (optically thin) boundary conditions discussed above, 
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the key question is how the outer disk solutions couple to the inner disk, where equation (10) is 
valid. In Appendix E, we show that at i? c i rc , the disk must satisfy a jump condition^], given by 

AW 3 

-Rcirc t -r 7- = — (25) 

W r] 

where AW' = W' + — W'_ is the jump in W' at -R c irc {i.e., the difference in the values of W at radii 
infinitesimally larger and smaller than i? c irc)- This result was first derived by J.E. Pringle (1997, 
private communication). Note that equation (25) indicates that W' + is larger than W'_; thus if 
W' + < 0, then W'_, the gradient of W just interior to i? C i rc , must be steeper (more negative). 

The outer disk equation can be solved numerically to find the appropriate boundary condition 
at -Rciro and then determine the matching inner disk solution, as in §4.1. However, we can 
show more directly that the solutions will always be close to the zero-crossing solutions. It is 
straightforward to find the gradient of the asymptotic solution for the outer disk, in the usual 
WKB approximation (e.g., Pender &: Orszag 1978 ) that W = e s (see also Appendix A): 



W' (5 + i\x) i 
X W = 2 2 



2 2 nr ■ 2 J 10 . A 2 2-25 f~,~ 3 

\x - 25 - ix \ h 4x x chc [ a + uj x° 



1/2 

(26) 



where x = x / x chc > 1- Comparison with numerical solutions of the outer disk equation shows 
that equation (26) (with choice of minus sign) is generally accurate to 10%, with the worst error 
being 15 — 20% (the latter being the case for some of the prograde disk modes, which tend to 
have smaller values of x c j rc than the retrograde modes); equation (26) always underestimates the 
magnitude of the gradient. (Note that equation [26] does not predict W' — > as x — > a?out; since 
the assumption that S" <C {S') 2 breaks down as x approaches x out .) Since we are chiefly interested 
in equation (26) for deriving the jump condition at the circularization radius, we take the limit 
x — ► x c i rc , so that x ~ * 1 ; 



W = (5 + ix chc ) 
X W 2 



x 2 

""are 



25 - i { 10x circ + Ax 4 c r 25 (or + cD 4 rc ) }] V2 . (27) 



In Figure 13 we plot the value of the gradient in the disk tilt, xW'/W, just outside the 
circularization radius, for all of the X-ray binary warped disk models presented by Maloney 
& Begelman (1997). The solid portions of the curves indicate retrograde modes. In all cases 
the gradient is negative and generally large, especially for the retrograde modes; to satisfy the 
jump condition, the gradient just inside -R c irc must be even steeper (note that xAW'/W = 6/77). 
Because the gradient is always negative and steep where the inner disk solution patches on to the 
outer disk solution, the former must have a shape close to that of the zero-crossing mode with a 



1 This jump condition is derived by assuming that mass and angular momentum are added via the accretion stream 
at a fixed radius. It assumes no viscosity discontinuities due to the stream, and hence is itself an idealization. On the 
other hand, assuming the warp goes to zero at -R c h- C is in some sense the opposite limit, as it implicitly posits that the 
effect of viscosity is to pin the disk to the orbital plane at the circularization radius. That these two extreme views 
yield comparable results is an indication of the insensitivity of our results to the exact choice of boundary conditions 
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zero at i? c i rc . Moreover, the tilt of the outer disk solution always goes rapidly to zero for R > i? c irc- 
For example, Figure 14 plots the tilt of the outer disk solutions as a function of radius, for the 
most rapidly precessing, steady-state (<7j = 0) solutions, for 5 = 0.75 and 1.25. Even for the most 
slowly-declining solution, the amplitude of the tilt at x ou t is less than 1% of the value at x c i rc . 

Changing the outer boundary condition does not allow a range of solutions to exist, unlike the 
optically-thin-boundary condition discussed above. The relevant physical question then becomes, 
what is the nature of the unstable mode for a given value of x c i rc and uj 1 In Figure 15 we show 
the precession and growth rates for the unstable modes as a function of u , for a fixed value of 
2xirc = 30. The upper half of the figure shows the precession rate and the lower half shows the 
growth rate. As expected, at small values of the quadrupole torque parameter the unstable mode 
is prograde. With increasing uj the unstable mode switches to retrograde precession, and the 
growth rate starts to decline, eventually going to zero. This is easy to understand physically: for 
a fixed value of x c i rc , the effect of increasing uj is to make the external torque more important at 
x c i rc , until eventually this dominates over the radiation torque in controlling the precession, and 
forces the warping mode to become retrograde. At very large values of the external torque, the 
viscous torques produced by the driven differential precession become too strong compared to the 
radiation pressure torques, and the instability is switched off; thus the growth rate goes to zero as 
uj approaches this critical value. 

There is actually a second branch of retrograde solutions displayed in Figure 15 (the dashed 
curve); these have comparable growth rates to the modes just discussed but precession rates about 
an order of magnitude larger. This second branch consists of modes for which the gradient in the 
tilt is positive (and generally large) at the circularization radius, unlike the solutions displayed in 
Figures 13 and 14. For these solutions the warp actually peaks in the outer disk at R > R c i rc 
and then declines rapidly towards zero. It is not clear that these solutions will ever occur in real 
astrophysical disks, as our treatment of the outer disk (.Rcirc < R < Rout) is much more of an 
idealization than that of the inner disk. 

5. Discussion 

Earlier work on the radiation-driven warping instability discovered by Pringle (P96; MBP) 
considered only the isothermal, 5 = 3/2 case. In this paper we have considered more general 
power-law disk density distributions, from the isothermal disk to 5 = —3/2, corresponding to a 
radiation-pressure supported disk; this spans the range that is likely to be relevant to astrophysical 
disks. Although the shapes of the eigenfunctions do change with decreasing 6, the most important 
features of the instability are generic. Most importantly, the instability exists over the entire range 
of surface density index that we have considered, and the critical radius above which disks are 
unstable to radiation-driven warping changes only by a factor of ~ 6 from 5 = 3/2 to 5 = —3/2. 
Similarly, the growth and precession rates (in dimensional units) do not depend strongly on S (see 
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the discussion after equation [18] and below). Evaluating equation (15) for the critical radius, 
R cr = f 5.9 x 10 8 - 3.5 x W 9 ) ( —\ cm, 



(5.9 x 10 



aKv ml^ ( V W M 



3.5 xlO 17 -*- cm, (28) 

* 1 e .i / UO 8 M ,' 



where the range in numerical values is for 5 = 3/2 to 5 = —3/2 and e = O.leo.i- The only 
warping modes with zeros at R < R cr are damped, so that disks that are smaller than R cr are 
stable against warping. In consequence of the e -2 scaling of R cr , accretion disks in systems with 
very low radiative efficiency will not be unstable to radiation-driven warping unless they are 
implausibly large. For this reason, this mechanism cannot provide an explanation for the warp 
in the thin maser disk of NGC 4258 (e.g., [Miyoshi et al. 1995j |Herrnstein et al. 1997|) if the 



inner disk is advection-dominated with e ~ 10 -3 (Lasota et al. 1996; see the discussion in MBP), 
since the maser disk would be far too small for instability in this case. This also indicates that 
radiation-driven warping generally will not be important in cataclysmic variables or protostellar 
disks dominated by accretion-powered luminosity, since the radiative efficiency is limited to small 



values as the stellar surfaces are at i?* 3> R s (but see Armitage fc Pringle 1997 for a discussion of 



the possible action of the instability in the protostellar case). 

To evaluate the typical precession timescales, we need to evaluate the viscous inflow timescale 
at R cr . Letting v\ = ac s H, where c s is the isothermal sound speed and H is the scale height, we 
can write the viscous timescale as 

t vl sc~\^-\H/R)- 2 (29) 

where Vj, is the rotational velocity (assumed to be Keplerian) and H/R is evaluated at the radius 
in question. Since t v isc oc i? 3/2 , and R CT /R 

o — x 



2 

cr ) 



W*,*) ~ I f-V — *4 /2 (H/R)- 2 (30) 



3 \ e / ac 

where H/R is now evaluated at R cr . Taking the precession timescale i prec = 2n/(j T , where a r is 
given by equation (18), and evaluating the constants, we find 



„, y 2 M/M Q (H/R\ 

i prec ~ 12 -3 — — days (31) 

eo.i «o.i VO.01/^ 



with only weak dependence on 8: the numerical coefficient only varies by a factor of two over the 
whole range of 5. Thus the precession timescales for X-ray binary systems (the only systems in 
which precession can actually be observed) are expected to be of the order of weeks to months. 

This is of course the precession timescale for the steady-state modes from linear theory. As 
discussed in §3, real disks will ordinarily be unwarped beyond some maximum radius, either the 
physical edge of the disk or where the disk becomes optically thin. This outer boundary, which will 
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not in general correspond to the critical radius, will determine the warp growth rate. We expect 
that the warp will eventually saturate at some amplitude (but see Pringle 1997). Assuming that 
the disk does reach a steady state, what will the precession rate be? There is reason to suspect 
it may not be very different from the linear theory result. Figure 6 shows that, except for growth 
rates very close to the maximum, the real part of the eigenvalue ay, i.e., the precession rate, is 
nearly independent of the growth rate. In the isothermal case, in fact, ay is independent of a^. 
This suggests that, however different modes may couple in reaching the final state, the precession 
rate will be similar to the linear steady-state result. 

Implicit throughout this paper has been the assumption that the disks are optically thick 
to both absorption and re-emission, so that they are subject to the radiation-driven warping 
instability. This requirement imposes a minimum mass accretion rate that must be exceeded for 
the disk to be optically thick. In Appendix D, we derive this critical mass accretion rate for three 
different possible sources of opacity in astrophysical disks (electron scattering, dust absorption, 
and Kramer's opacity) and show that it does not in general place any significant limitations on 
occurrence of the instability. 

As discussed in §3.2, there is one very important systematic change in the nature of the 
instability with 5. The difference in the behavior of the growing modes for S > 1 and 5 < 1 is of 
fundamental importance for the evolution of disks warped by radiation pressure. For 6 > 1, the 
fast-growing modes all have their maximum warp (i.e., tilt (5) close to the outer edge of the disk, 
and the amplitude /3 max is much greater than (3 , the tilt at the origin. This immediately implies 
that the warp must reach the nonlinear regime when the tilt at small radius is negligible. In this 
case the evolution of the disk at radii interior to the warp maximum is almost certainly driven by 
the nonlinear evolution of the outer warp (e.g., |Pringle 1997] ), so that nonlinear effects must be 
important if the warp extends to the disk inner edge. 

For 5 < 1, the behavior is qualitatively different, as f3 max /f3 is always of order unity. In 
this regime, nonlinearity will be important only if the warp has grown out of the linear regime 
at the origin. Furthermore, because the shapes of the growing warps in these two regimes are so 
dissimilar, the effects of shadowing of the central source by the warping of the disk will be very 
different. These distinctions are liable to be crucial for X-ray binaries such as SS 433 and Her X-l, 
which show evidence for a global precessing warp. 

One final point regarding X-ray binary systems must be mentioned. In one of the best-studied 
systems, Her X-l, the direction of precession of the warp is inferred to be retrograde with respect 



to the direction of rotation (e.g., Gerend fc Boynton 1976 ) and this has also been suggested for 



SS 433 (Leibowitz 1984; Brinkmann, Kawai & Matsuoka 1989). As shown in Appendix B, in the 



absence of external torques the direction of precession of the warp must be prograde. However, 
the qualification on this statement is extremely important: as pointed out in §4.2, and discussed 
in detail by Maloney & Begelman (1997b), including the quadrupole torque from a companion 
star allows retrograde as well as prograde solutions to exist. 
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The zero-crossing outer boundary condition that we have imposed will not be strictly correct 
in real astrophysical disks. However, as discussed in §4, the solutions that obey the likely realistic 
outer boundary conditions - the optically thin outer boundary for accretion disks in active galactic 
nuclei, and a flat outer boundary for disks in X-ray binaries - are in all important respects similar 
to the zero-crossing solutions. 

Radiation-driven warping and precession offers a robust mechanism for producing tilted, 
precessing accretion disks, in accreting binary systems such as Her X-l and SS 433, and in active 
galactic nuclei. Because radiation-driven warping is an inherently global mechanism, it avoids the 
difficulties inherent in other proposed mechanisms for producing warping and precession, e.g., 
communicating a single precession frequency through a fluid, differentially-rotating disk. This 
mechanism can thus explain the simultaneous precession of inner disks (as evidenced by the jets 
of SS 433 and the pulse profile variations of Her X-l) and outer disks (as required to match the 
periodicities in X-ray flux and disk emission in these same objects). 

A full understanding of the nature of the radiation-driven warping instability will require 
nonlinear simulations of the type presented in Pringle (1997), which will not only allow for 
inclusion of the nonlinear terms but also inherently nonlinear effects such as shadowing. This will 
be the subject of future work. 
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to Jim Pringle for his insightful comments regarding the choice of outer boundary conditions. We 
would also like to thank the referee for helpful comments on the paper. PRM was supported by 
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support from NSF grant AST-9529175. The research of PRM and MCB was supported by NASA 
grant NAG5-4061 under the Astrophysical Theory Program. MAN was supported by the NASA 
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A. Asymptotic Solutions to the Twist Equation 

In this appendix we examine the behavior of the twist equation at large x, and show that, 
unlike the 5 = 3/2 case studied by MBP, the solutions for W always diverge as x — > oo. We also 
show from the asymptotic solutions that the eigenfunctions for the case 5 = 1 always diverge as 
the growth rate — <jj — > 0.25, independent of the real eigenvalue ay, as is seen in the numerical 
solutions (§3.2). 

We write the twist equation (10) as 

xW" + (2 - ix)W' - idx^W = (Al) 
where primes denote derivatives with respect to x and // = 3 — 25 > 0for5< 3/2; recall that x is 
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real and a is complex. For any value of /i, this equation has an irregular singular point at infinity. 
Assuming that W = e s (e.g., Bender Orszag 1978[ ), we obtain 

xS" + x (S') 2 + (2 - ix)S' - iax^ = 0. (A2) 

We further assume that S" < (S'y. As x — > oo, the 2S' term can be neglected compared to 
—ixS'; however, we retain it for reasons that will become obvious below. Equation (A2) is then 

(S') 2 + - i\ S' - iax K = (A3) 



where k = [a — 1 > — 1. The solution to this equation is 

\ 1/21 

(A4) 



2 



x 



This equation has qualitatively different behavior at large x for k < and k > 0, corresponding to 
5 > 1 and <5 < 1, respectively. (This difference in behavior is of course also seen in the numerical 
solutions, as discussed in §3.) If k < 0, then all the x— terms in the square root are small compared 
to unity, and expansion of the square root gives the two solutions 

S f ~i-~ + ax K , -ax K ; -1 < k < 0. (A5) 

x 

(Note that the validity of the assumption that S" -C {S') 2 requires that ax K 3> k/x; this obviously 
breaks down as k — > — 1.) Integrating and exponentiating then gives the two solutions for W 

W + ~ e ix x- 2 e^^ )x \ W.-e-^ 1 ^; < n < 1 (A6) 
W + ~ e"x^ 2 , W_ ~ x" 5 ; n = 0. (A7) 

The general solution will be a linear combination of W+ and Equation (A7) is just the 

asymptotic behavior of the Kummer functions M(a,b,x) (cf. MBP; Abramowitz & Stegun 1964, 
eq. [13.5.1]), with a = a and 6 = 2. Equation (A7) shows that for /i = we can impose the 
condition W — > as x — > oo provided that < a r < 2. However, equation (A6) shows that it is 
not possible to impose this condition for < fi < 1 if a r ^ 0: one solution for W always diverges 
as x — ► oo. Thus, unlike the degenerate 5 = 3/2 solutions, W always diverges as x — > cx>, for any 
value of & r ■ 

For k = 0, corresponding to 5 = 1, the a— term under the square root in equation (A4) is 
independent of x. Since \&ia\ is not necessarily small compared to unity, we expand the square 
root as 

\ _ m + fiV /2 ~ U _ 1_) (A8) 



x j \ <jjx 

with u = l — Ma. We thus obtain 



S±~^-\nx±uj 1 / 2 (ix-1-lnx) . (A9) 
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The first term inside the parentheses on the righthand side of (A9) will cause W to diverge as 
x — ► oo, as exp(Im(u; 1//2 )x). Since Re(u; 1 / 2 ) — > as <7j — > —0.25, the second term will cause VF to 
diverge as <7j — > —0.25, independent of <5y (see also Appendix C). This divergence is reflected in 
the behavior of the zero-crossing solutions; see §3.2 and Figures 5a and 5b. 

Finally, consider the case k > 0. In this case, at large x the x K term under the square root of 
equation (A4) dominates, and so 

5 ± ~|-lnx±^^x A (A10) 

where A = 1 + k/2 = (/x + l)/2; A > 1 for 5 < 1. For any cr/ 0, the third term causes W to 
diverge exponentially as x — > oo, as for the other solutions for (5 < 3/2. 



B. Proof of the Distinctness of the Eigenvalues 

We separate the twist equation into real and imaginary parts: 

x/3" + xM (1 - 7) + 2/3' = -x^[5di (Bl) 
X1 " + -J (2Y - 1) + 2Y = s"er r (B2) 

where // = 3 — 2(5 and primes denote derivatives with respect to x. Equation (B2) shows that a 
necessary condition for an eigenvalue to occur is 27' — 1 = (i.e., 7' — > 1/2 as (3 — > 0). Multiplying 
equation (B2) by 2x/3 2 and rearranging, we obtain 



x 2 /? 2 (2 7 ' - 1) = 2xf3 2 (x^a r - 1) . (B3) 



Integration then gives 



x 2 /3 2 (2 7 ' - 1) = 2 I* /3 2 x (x^d r - 1) dx . (B4) 
•/ 

Consider first the case \i > (5 < 3/2). Changing variables to w = ayx^, we rewrite equation (B4) 
as 

The integral on the righthand side of (B5) is monotonic; hence, there can be at most one zero for 
a given a r . Note also that <5y must be greater than zero for there to be an eigenvalue, i.e., the 
precession of the warp must be prograde (in the same direction as the rotation of the disk) . 

The change of variables to u is not valid for /x = (<5 = 3/2). In this case equation (B3) is 
simply 

x 2 /3 2 (27 / -l)]' = 2x/3 2 (<7 r -l) . (B6) 
If ay = 1, the righthand side is identically zero, and so 

(B7) 
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where C is a constant. The boundary condition 7' — ► ay/2 as x — > for 5 = 3/2 (c/. equation [12] 
in the text) requires that C = 0, and thus for ay = 1, 7' = 1/2 for all x, as can be seen from the 
Kummer M functions with a = ay = 1 (MBP). 

If ay 7^ 1, we integrate equation (B6) to obtain 

( 2 V - !) = ( ? r " 1) r /? 2 ^ . (B8) 

X p Jo 

This equation shows that there are no zeros for ay ^ 1, as the integral on the righthand side is 
positive definite. Instead we simply have 7' > 1/2, ay > 1; 7' < 1/2, ay < 1. 



C. Solution for the case 5 = 1 



The twist equation (10) can be re- written as 



d_ 

dx 



x 2 e- ix 



dW 

dx 



+ x 2 e~ ix 



ix 2 ~ 25 aW 



which for 5 = 1 simplifies to 

Define the new function <f) by 
so that 



d_ 

dx 



x 2 e~ lx 



+ x 2 e~ ix [iaW] = 0. 



d, x 
W = x- x e ixl2 



^ = X-V X ' 2 ^ + 

OX ox 



Then the twist equation (C2) can be re-written as 



2 



dx 2 



+ 



1 ._\ i 
- - ic + - 
4 x 



2~ X 







(CI) 

(C2) 

(C3) 
(C4) 



(C5) 



If we now define b 2 = —(1 — Aicr) 1 and a new radial variable y = x/b, equation (C5) becomes 

d 2 



dy< 



+ 



1 ib 

■7 + — 
4 y 



. 



This is Whittaker's equation (e.g., Abramowitz & Stegun 1964, equation [13.1.31]), 



d 2 u 



dz 



+ 



1 K (l ~ 



A+- + 

4 2 



to = 



where the Whittaker functions w = Af K „ are related to the Kummer functions M(a, b, z) by 

M K) p = e-iz^M Q + M -k,1 + 2 M ,* 
where in our case // = 1/2, k = i6. In terms of W and x, this solution is 

f 1— ib) / X 

W(x) = e — ~ X M (l — ib,2,— 



(C6) 



(C7) 



(C8) 



(C9) 
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D. Critical Accretion Rate 

For an accretion- fueled source, neither the growth/precession rates nor the critical radius 
for instability explicitly depends on the mass accretion rate, but only on the radiative efficiency 
e = L/Mc 2 . However, the requirement that the disk be optically thick to absorption and 
re-emission of the incident flux does impose a minimum value for the surface density £ oc M/a in 
the usual a— viscosity formalism. 

We define a critical mass accretion rate M cr such that the disk is optically thick at the critical 
radius for instability, i.e., the surface density £ = S cr there. We consider three possible sources of 
opacity: 

(1) Dust opacity. This will dominate if the gas is molecular or atomic with T £ 10 4 K. (The 
dust will be considerably cooler than the gas.) In this case the Rosseland mean opacity 
coefficient kr = 0.7, 2.7, 5.0, and 7.5 cm 2 gm -1 for dust temperatures = 50, 100, 300, 
and 600 K, respectively ( Pollack et al. 1994 ; this assumes Solar neighborhood abundances 



and depletion patterns). The corresponding critical surface densities are 

S cr (dust) ~ 4.3 x lO 23 ^ 1 cm" 2 ~ 6 x 10 22 - 6 x 10 23 cm" 2 (Dl) 
for T d = 50 - 600 K. 

(2) Electron scattering. Ignoring Klein-Nishina corrections, in this case kr ~ 0.35 cm 2 gm , 
and so 

S cr ~ 1.2 x 10 24 cm" 2 (D2) 

independent of n or T. 



(3) Kramer's (bound-free and free-free) opacity for ionized gas. From Schwarzschild (1958) the 
Rosseland mean opacity in this case is approximately kr ~ 3 x 10 23 pT" 7 / 2 cm 2 gm -1 for 
electron temperatures T ^ 10 6 -fT, which gives a critical density 

S cr « 6 x 10 23 n#r 7/2 cm" 2 . (D3) 

The numerical coefficient is for T ~ 10 6 K; it slowly decreases with increasing T. 

From the definition of R Q (cf. equation [8]) we have 

4 (D4) 



i? cr > 



R a 

and so the requirement that the disk be optically thick at the critical radius then becomes 

5^cr = ^o^cr ^ (D5) 
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which can be written as a constraint on S G = £(i? ): 

S D = £ cr x 2 r 5 . (D6) 

The surface density T, = M /3ttv , which in the a— viscosity prescription is M /3ira(c s H) , where 
the isothermal sound speed c s and the disk scale height H are evaluated at R Q . For a thin 
disk with negligible self-gravity, H ~ Rc s /V OT b. From the definition of R Q , V 0T b(R ) = {e/rfjc. 
To adequate accuracy we can take the sound speed to be c s ~ 13T^ 2 km s -1 , where the gas 
temperature T = 10 4 T4 K, and so 



,M 



3 



S = 3.8 x 1(T 3 — {TiRsy 1 ^-J . (D7) 
where R s is in units of cm. If we write M in terms of the Eddington accretion rate 

M e e^ = 4.7x W U R S gm s" 1 , (D8) 
& 

then 

£ = 7.7 x 10 32 — T 4 _1 (- Y cm" 2 (D9) 
a \r}J 

where m = M /Me- Equating this to S cr x 2 J?, we finally obtain an expression for the critical mass 
accretion rate 

- = 1-3 x 10- 10 -|^T 4 ( V -)\^. (DID) 
a W Zii cm z \e J 

Equation (D10) is simply applicable for the case of dust or electron scattering opacity. However, 
for Kramer's opacity, we need to know the density as well as the column density. Solving for the 
density at R Q and proceeding as before, we obtain a slightly different equation for the critical 
accretion rate: 

9/2 

a " " \ ej 

where Tq = T/10 6 K, a more realistic value for the Kramer's opacity case. 

Consider first equation (D10), for the case of dust and Thompson opacity: 

(1) The dust temperature is likely to be Td ~ a few hundred K when this is the dominant source 
of opacity, so that S cr ps 10 23 cm -2 , while T 4 ~ 0.1 - 1. For 5 = 3/2 this is 

3 

— = 3.2 x KT 5 T 4 4- (D12) 

where x CT = 2ir and we have normalized e = O.leo.i- For 5 = 1, the numerical coefficient 
is reduced by a factor of 5. This does not impose a serious constraint on disk stability; of 
course, ih cr oc x 26 , so for 5 > 0, opacity effects become more stringent if the disk is to be 
optically thick at R R cr (x » x cr ). 



/ \ 9/2 

- = 7.6xl0- 11 r 3 i?y 2 (^) x s cr (Dll) 
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(2) For electron scattering S cr ~ 1.2 x 10 24 cm 2 , while T }t 10 6 K, so 



m - 1.6 x 10-*T 6 4-x 2 J. (D13) 



=0.1 



This requirement is most stringent for 5 > 0, as x„ ranges from one to (27r) 3 ~ 248 for 
< 5 < 3/2, although an isothermal disk is unlikely to be a good assumption in the 
electron-scattering regime. 

Now consider equation (Dll), for Kramer's opacity. 

(3) In this case Tq jt 1, and substituting for R s in equation (Dll) 

= 1.3X10-T|( ) ^ (D14) 
x u/ e o.i 

where M is the mass of the central object. For active galactic nuclei this can be a significant 
constraint if 6 > 0, since M ~ 10 6 — 10 9 Mq. However, Tg <C 1 for accretion disks around 
AGN at the radii of interest for warping, so this is probably irrelevant. 



E. Disk Outer Boundary Jump Condition in X-Ray Binaries 



We start with the basic conservation equations from Pringle (1992) - his equations (2.1) 
and (2.2). We assume that matter is injected into the disk at i? c irc, and that Vr is zero for 
-Rtirc < R < -Rout- Integrating the mass and angular momentum conservation equations across 
-Rtirc gives the jump conditions 



Rare 

T/ n '\l 1 V 91 



M 



Rr 



(El) 



(E2) 



where 1 is the unit tilt vector giving the direction of angular momentum at R, as before, and 1 Q is 
the unit vector normal to the equatorial plane of the binary: 



lo = (0,0,1) 



(E3) 



in our usual Cartesian coordinate system. Plus and minus signs denote quantities evaluated 
at radii just outside and inside of i? c iro respectively. We assume that the disk is continuous, 
so that 1+ = 1_ = 1, but that the derivative dl/dR may be discontinuous across the boundary. 
Furthermore, since there may also be a discontinuity in the surface density, we assume that 



(Sz/i), = (Sz/i)_ 



5d 
S 



(£i/i)_y. 



(E4) 
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Substituting from equation (El) into (E2) and rearranging, we have 



V 



d\\ / 31 \ 

dRJ, ~ [dRJ 



M 

(1 - lo) • (E5) 



-Rcirc (S^l) 



Given (3'_, and (S^i) , equation (E5) gives the three conditions necessary to determine y, f3' + , 
and r y' + . For example, to find y we take the dot product of equation (E5) with 1: since 



dRJ + \dRj_ 
we have 

fe-^-Tdkr* 1 -™ 3 ' <E7) 

which gives the useful result that in the small angle (/3 <C 1) limit, y = 1 + O (/3 2 ), so that the 
surface density is continuous in this regime. Similarly, in this limit the Z-component of equation 
(E5) is also O (/3 2 ), and so we neglect it. In terms of W = /3e n , equation (E5) in the linear (small 
(3) regime then becomes 

AW 2M 



Wa TC 7]R ciTC (Sl/i) 

Since in the small angle limit we also have 



(E8) 



(EV R )_ = - (E^)_ - = -i-^= (E9) 



Reive \ ^ / 2 Rcirt 



(assuming a Keplerian rotation curve), we finally obtain 

AW' 3 

itWc^ = - (E10) 

Wcirc V 

which is the desired jump condition. Equation (E10) was first derived by J.E. Pringle (1997, 
private communication) . 
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Fig. 1. — The location of the zero, x Q , for the first ten zero-crossing eigenfunctions, for the steady- 
state (crj = 0) solutions of the twist equation. The locations of the zeros are plotted as a function 
of the surface density power-law index 5, for —3/2 < 5 < 3/2. The steps in d are 0.05. The merging 
of eigenvalues that occurs for the higher-order eigenfunctions at S ^ 1 is marked. 

Fig. 2. — The real eigenvalues for the same steady-state eigenfunctions as in Figure 1. Plotted are 
the real eigenvalues ay for the first ten zero-crossing eigenfunctions. The case 6 = 3/2 is obviously 
degenerate (ay = 1 for all order eigenfunctions). The first-order eigenfunction has the largest 
eigenvalue for all other values of 5. The behavior of the higher-order (more distant) zeros becomes 
very complex for smaller values of S, with merging of the eigenvalues, as seen in Figure 1. 

Fig. 3. — The magnitude of the tilt f3 as a function of x, for the steady-state eigenfunctions. Plotted 
is /3(x) for (left to right) 6 = 1.5 to —1.5, in steps of 0.25 in 5. The tilt at the origin has been set 
to unity. 

Fig. 4. — The location of the zero for the (growing) time-dependent eigenfunctions. Plotted is x for 
the first-order growing modes, as a function of the growth rate — oV From left to right, the curves 
are for 5 = —3/2 to S = 1.45, in steps of 0.05 in S. The upper envelope delineates the maximum 
growth rates for the first-order modes. For S = 1, x — > oo as — di — > 0.25; see Appendices B and 
C. The curve for this case has been plotted to — d\ = 0.249. 

Fig. 5. — (a) The location of the zero as a function of growth rate — di, for the first three order 
eigenfunctions for 5=1. The zeros merge as — di — ► 0.25 (see Appendices B and C). (b) The real 
part of the eigenvalue, function of growth rate — di, for the three 5=1 eigenfunctions of 

Figure 5a. The real part of the eigenvalue approaches zero as — di — > 0.25; see Appendices B and 
C. 

Fig. 6. — The real part of the eigenvalue, function of growth rate — di, for the first-order 

growing modes shown in Figure 4. For <5 = 1, ay — > as — di — > 0; for 5 = 3/2, a r is independent 
of di. 

Fig. 7. — The tilt (3 as a function of radius variable x for the fastest-growing modes. The maximum 
value of the tilt has been normalized to unity. Plotted are (from left to right) (3{x) for 5 = 1.45, 
1.35, 1.25, 1.15, and 1.05. 

Fig. 8.— 0(x), as in Figure 7, for (from left to right) S = 0.55, 0.65, 0.75, 0.85, and 0.95. As in 
Figure 7, the maximum value of (5 has been normalized to unity. 

Fig. 9. — The real portion of the eigenvalue ay versus the imaginary part of the eigenvalue aj, 
for warping modes that satisfy the optically thin outer boundary condition discussed in §4.1. The 
optically thin radius has been set equal to the radius of the first-order zero for the steady-state 
solutions for (from left to right) 5 = 1.50, 1.25, and 0.75. The dotted portions of the curves have 
di > 0, and are therefore damped. 
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Fig. 10. — The tilt (3 as a function of radius variable x, for the steady-state zero-crossing modes 
(dotted lines) and the fastest-growing optically-thin-boundary modes from Figure 9. The tilt has 
been normalized to a maximum of unity in all cases. 

Fig. 11. — As in Figure 9, but with i? t hin set equal to the zero for rapidly- growing zero-crossing 
modes rather than steady-state zero-crossing modes. 

Fig. 12. — As in Figure 10, but with i? t hin set equal to the zero for rapidly-growing zero-crossing 
modes rather than steady-state zero-crossing modes. 

Fig. 13. — The gradient in the disk tilt, xW'/W, just outside the circularization radius, for all of 
the prograde and retrograde warp models presented in Maloney & Begelman (1997b). This has 
been derived using the asymptotic solution to the outer disk equation, given by equation (27). The 
solid lines indicate the solution regimes where the precession direction is retrograde; the prograde 
solutions are: 6 = 1.25, <7j = (dotted line); 5 = 0.75, &i = (short-dashed line); 5 = 0.75, 
maximum — U{ (long-dashed line); 5 = 1.25, maximum — U{ (dot-dashed line). Note that in all cases 
the gradient is large and negative. 

Fig. 14. — The tilt (3 of the outer disk solutions, for the most rapidly-precessing, steady-state warp 
modes found by Maloney & Begelman (1997b), for 5 = 0.75 and 1.25. The retrograde solutions are 
shown by the solid (S = 1.25) and dotted (5 = 0.75) lines, while the prograde solutions are shown 
by the short-dashed (8 = 1.25) and long-dashed (5 = 0.75) lines. 

Fig. 15. — The steady-state and growing modes for a fixed value of the circularization radius, 
^circ = 30. The disk truncation radius x out is assumed to be \/3 times larger. The precession 
rate a r (top panel) and growth rate di (bottom panel) are plotted as functions of the companion 
torque parameter cj . The dotted curve shows the prograde modes and the solid curve shows the 
retrograde modes. The dashed line is a second branch of retrograde solutions; it is not clear that 
these will exist in real accretion disks (see discussion in §4.2). 
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